library(survey)
library(foreign)
library(ggplot2)
library(cowplot)
library(dplyr)

gss.df <- read.dta("~/Downloads/GSS7218_R2.DTA")

#### CLEANING ####

## run cleaning from GSS_figures.R

## Confidence in other institutions
gss.df$conscotus <- ifelse(as.numeric(gss.df$conjudge) == 1, 1, 0)
gss.df$conscotus[gss.df$conjudge == "NA"] <- NA

gss.df$concongress <- ifelse(as.numeric(gss.df$conlegis) == 1, 1, 0)
gss.df$concongress[gss.df$conlegis == "NA"] <- NA

gss.df$conmedia <- ifelse(as.numeric(gss.df$conpress) == 1, 1, 0)
gss.df$conmedia[gss.df$conpress == "NA"] <- NA

gss.df$conbank <- ifelse(as.numeric(gss.df$confinan) == 1, 1, 0)
gss.df$conbank[gss.df$confinan == "NA"] <- NA

gss.df$conpres <- ifelse(as.numeric(gss.df$confed) == 1, 1, 0)
gss.df$conpres[gss.df$confed == "NA"] <- NA

gss.df$conorglabor <- ifelse(as.numeric(gss.df$conlabor) == 1, 1, 0)
gss.df$conorglabor[gss.df$conlabor == "NA"] <- NA

gss.df$conbigbus <- ifelse(as.numeric(gss.df$conbus) == 1, 1, 0)
gss.df$conbigbus[gss.df$conbus == "NA"] <- NA

gss.df$conrelig <- ifelse(as.numeric(gss.df$conclerg) == 1, 1, 0)
gss.df$conrelig[gss.df$conclerg == "NA"] <- NA

gss.df$conedu <- ifelse(as.numeric(gss.df$coneduc) == 1, 1, 0)
gss.df$conedu[gss.df$coneduc == "NA"] <- NA

gss.df$conmedicine <- ifelse(as.numeric(gss.df$conmedic) == 1, 1, 0)
gss.df$conmedicine[gss.df$conmedic == "NA"] <- NA

gss.df$contele <- ifelse(as.numeric(gss.df$contv) == 1, 1, 0)
gss.df$contele[gss.df$contv == "NA"] <- NA

gss.df$conscience <- ifelse(as.numeric(gss.df$consci) == 1, 1, 0)
gss.df$conscience[gss.df$consci == "NA"] <- NA

## Other covariates for Figure 6.3 (logit model)

gss.df$black2 <- ifelse(gss.df$race == "black", 1, 0)
gss.df$waryear <- ifelse(gss.df$year == 1991 | gss.df$year > 2000, 1, 0)
gss.df$confidence <- (gss.df$conscotus + gss.df$concongress + gss.df$conmedia +
                        gss.df$conpres + gss.df$conedu + gss.df$conmedicine +
                        gss.df$conscience)/7
gss.df$dem2 <- ifelse(as.numeric(gss.df$partyid) < 4, 1, 0)
gss.df$rep2 <- ifelse(as.numeric(gss.df$partyid) > 4 & as.numeric(gss.df$partyid) < 8, 1, 0)

## Create weighted survey design object
gss_sd <-
  svydesign(
    id = ~ 1,
    weights = ~ wtssall,
    data = gss.df
  )

#### FIGURE 6.1 - CONFIDENCE OVER TIME ####

conf1 <- svyby(~ con_mil + conscotus + concongress + conpres +
                 conmedia + conscience + conedu + conmedicine, 
               ~ year, gss_sd, svymean, na.rm = TRUE)
conf1 <- conf1[-c(1,12), -c(10:17)]
for (i in 2:9) {
  conf1[,i] <- conf1[,i] * 100
}

insts <- c("Military"="solid", "SCOTUS"="dotdash","Congress"="longdash",
           "POTUS"="twodash", "Science"="dotted", "Education"="dashed",
           "Media"="12345678", "Medicine"="4C88C488")
c1_plot <- ggplot(conf1) +
  geom_line(aes(x = year, y = con_mil, linetype = "Military")) +
  geom_line(aes(x = year, y = conscotus, linetype = "SCOTUS")) +
  geom_line(aes(x = year, y = concongress, linetype = "Congress")) +
  geom_line(aes(x = year, y = conpres, linetype = "POTUS")) +
  geom_line(aes(x = year, y = conscience, linetype = "Science")) +
  geom_line(aes(x = year, y = conedu, linetype = "Education")) +
  geom_line(aes(x = year, y = conmedia, linetype = "Media")) +
  geom_line(aes(x = year, y = conmedicine, linetype = "Medicine")) +
  labs(x = "Year",
       y = "Public with Confidence in Institution (%)",
       linetype = "Institutions") +
  xlim(1970,2020) + ylim(0,65) +
  scale_linetype_manual(values = insts) + theme_bw()
c1_plot
ggsave("fig6.1.jpeg", plot = c1_plot, dpi = 500)

#### FIGURE 6.2 - CONFIDENCE OVER TIME, MEAN ####

conf2 <- conf1
conf2 <- conf2 %>%
  mutate(conf_means = rowMeans(select(., c(conscotus, concongress, conpres, conscience,
                                           conedu, conmedia, conmedicine))))
conf2 <- conf2[,-c(3:9)]

insts2 <- c("Military"="solid", "Other Institutions"="dashed")
c2_plot <- ggplot(conf2) +
  geom_line(aes(x = year, y = con_mil, linetype = "Military")) +
  geom_line(aes(x = year, y = conf_means, linetype = "Other Institutions")) +
  labs(x = "Year",
       y = "Public with Confidence in Institution (%)",
       linetype = "Institutions") +
  xlim(1970,2020) + ylim(20,65) +
  scale_linetype_manual(values = insts2) + theme_bw()
c2_plot
ggsave("fig6.2.jpeg", plot = c2_plot, dpi = 500)

#### FIGURE 6.3 - PREDICTED CONFIDENCE AGAINST OTHER INSTITUTIONS ####

formula63 <- con_mil ~ confidence + black2 + male + education + age + dem2 + rep2 + 
  waryear + year
fig3_mod <- svyglm(formula63, gss_sd, family = binomial(link = "logit"))
summary(fig3_mod)

pred_1974 <- data.frame(confidence = seq(0, 1, 0.1),
                        black2 = rep(0, 11),
                        male = rep(1, 11),
                        education = rep(mean(gss.df$education, na.rm = TRUE), 11),
                        age = rep(mean(gss.df$age, na.rm = TRUE), 11),
                        dem2 = rep(0, 11),
                        rep2 = rep(0, 11),
                        waryear = rep(0, 11),
                        year = rep(1974, 11))
preds74 <- as.data.frame(predict(fig3_mod, pred_1974, type = "response", se.fit = TRUE))

pred_2018 <- data.frame(confidence = seq(0, 1, 0.1),
                        black2 = rep(0, 11),
                        male = rep(1, 11),
                        education = rep(mean(gss.df$education, na.rm = TRUE), 11),
                        age = rep(mean(gss.df$age, na.rm = TRUE), 11),
                        dem2 = rep(0, 11),
                        rep2 = rep(0, 11),
                        waryear = rep(1, 11),
                        year = rep(2018, 11))
preds18 <- as.data.frame(predict(fig3_mod, pred_2018, type = "response", se.fit = TRUE))

hist74 <- hist(gss.df$confidence[gss.df$year == "1974"])
hist74$prop <- hist74$counts/sum(hist74$counts)*100
hist18 <- hist(gss.df$confidence[gss.df$year == "2018"])
hist18$prop <- hist18$counts/sum(hist18$counts)*100

conf3_74 <- data.frame(preds = preds74$response * 100,
                       ci_lo = (preds74$response - 1.96 * preds74$SE) * 100,
                       ci_hi = (preds74$response + 1.96 * preds74$SE) * 100,
                       confs = pred_1974$confidence,
                       hist = c(hist74$prop, NA))

conf3_18 <- data.frame(preds = preds18$response * 100,
                       ci_lo = (preds18$response - 1.96 * preds18$SE) * 100,
                       ci_hi = (preds18$response + 1.96 * preds18$SE) * 100,
                       confs = pred_2018$confidence,
                       hist = c(hist18$prop, NA))

conf3_plot74 <- ggplot(conf3_74) +
  geom_col(aes(x = confs, y = hist), fill = "gray80") +
  geom_line(aes(x = confs, y = ci_lo), linetype = "dashed") +
  geom_line(aes(x = confs, y = ci_hi), linetype = "dashed") +
  geom_line(aes(x = confs, y = preds)) + ggtitle("1974") +
  xlab("Mean Confidence in Other Institutions") + ylim(0,100) +
  ylab("Respondents with Confidence (%)") + theme_bw()
conf3_plot18 <- ggplot(conf3_18) +
  geom_col(aes(x = confs, y = hist), fill = "gray80") +
  geom_line(aes(x = confs, y = ci_lo), linetype = "dashed") +
  geom_line(aes(x = confs, y = ci_hi), linetype = "dashed") +
  geom_line(aes(x = confs, y = preds)) + ggtitle("2018") +
  xlab("Mean Confidence in Other Institutions") + ylim(0,100) +
  ylab("Respondents with Confidence (%)") + ylab("") + theme_bw()
conf3_plot <- plot_grid(conf3_plot74, conf3_plot18)
conf3_plot
ggsave("fig6.3.jpeg", plot = conf3_plot, dpi = 500)

#### FIGURE 6.4 - CONFIDENCE OVER TIME, BY PARTY ####

conf4 <- svyby(~ con_mil + conscotus + concongress + conpres + conscience + conedu +
                 conmedia + conmedicine, ~ year + dem, gss_sd, svymean, na.rm = TRUE)

conf4_rep <- conf4[c(1:32),]
conf4_rep <- conf4_rep[-c(1,12), -c(2,11:18)]
for (i in 2:9) {
  conf4_rep[,i] <- conf4_rep[,i] * 100
}
conf4_dem <- conf4[c(33:64),]
conf4_dem <- conf4_dem[-c(1,12), -c(2,11:18)]
for (i in 2:9) {
  conf4_dem[,i] <- conf4_dem[,i] * 100
}

conf4_plot_d <- ggplot(conf4_dem) +
  geom_line(aes(x = year, y = con_mil, linetype = "Military")) +
  geom_line(aes(x = year, y = conscotus, linetype = "SCOTUS")) +
  geom_line(aes(x = year, y = concongress, linetype = "Congress")) +
  geom_line(aes(x = year, y = conpres, linetype = "POTUS")) +
  geom_line(aes(x = year, y = conscience, linetype = "Science")) +
  geom_line(aes(x = year, y = conedu, linetype = "Education")) +
  geom_line(aes(x = year, y = conmedia, linetype = "Media")) +
  geom_line(aes(x = year, y = conmedicine, linetype = "Medicine")) +
  xlim(1970,2020) + ylim(0,80) + ggtitle("Democrats") + theme_bw() +
  xlab("Year") + ylab("Public with Confidence in Institutions (%)") +
  theme(legend.position = "none")
conf4_plot_r <- ggplot(conf4_rep) +
  geom_line(aes(x = year, y = con_mil, linetype = "Military")) +
  geom_line(aes(x = year, y = conscotus, linetype = "SCOTUS")) +
  geom_line(aes(x = year, y = concongress, linetype = "Congress")) +
  geom_line(aes(x = year, y = conpres, linetype = "POTUS")) +
  geom_line(aes(x = year, y = conscience, linetype = "Science")) +
  geom_line(aes(x = year, y = conedu, linetype = "Education")) +
  geom_line(aes(x = year, y = conmedia, linetype = "Media")) +
  geom_line(aes(x = year, y = conmedicine, linetype = "Medicine")) +
  xlim(1970,2020) + ylim(0,80) + ggtitle("Republicans") + theme_bw() +
  scale_linetype_manual(values = insts) +
  labs(x = "Year", y = "", linetype = "Institutions")
conf4_plot <- plot_grid(conf4_plot_d, conf4_plot_r,
                        rel_widths = c(1,1.45))
conf4_plot
ggsave("fig6.4.jpeg", plot = conf4_plot, dpi = 500)
